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An algorithm, complete with a specific FORTRAN implementation, is presented for the problem 
of determining whether or not two undirected graphs are isomorphic. The algorithm, centered upon 
the eigenvalues and eigenvectors of a modified adjacency matrix and techniques for decreasing the 
size of the automorphism group, is quite different from others (most of which are combinatorially 
based) and tends to work relatively very quickly on difficult test cases as well as on typical 
examples. Complexity estimates are given for many eventualities. 
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1. Introduction 

Two undirected graphs are said to be isomorphic if there exists a one-to-one correspondence 
between their nodes which preserves adjacency. A computationally difficult issue in the theory of graphs 
is the determination of isomorphism. For example the graphs in figure la are pairwise isomorphic while 
those in lb are not. 
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The ability to divide a set of graphs into maximal subsets whose elements are pairwise isomorphic 
has several immediate applications. If nodes are identified with atoms of a particular compound and if 
edges are identified with the existence of a chemical bond, then, for instance, such an ability is of use 
to the chemist in distinguishing among many theoretically produced compounds. Other applications 
include physics (e.g. "hearing" the shape of a drum [1] ) 1 , electronic circuit theory, linguistics etc. 

Much effort has been devoted to the development of an algorithm which will determine whether or 
not any two graphs are isomorphic in an amount of time that is bounded by a finite polynomial in the 
number of nodes of the graphs considered. At the moment, the authors are unaware of the existence of 
any such algorithm. Many approximate or heuristic procedures have been developed however. Many 
such algorithms check conditions combinatorially necessary for isomorphism, some also attempt to 
construct a permutation that would exhibit isomorphism. Corneil presents a review of many such tech- 
niques. All known exact methods grow exponentially in required time. Unfortunately, as pointed out in 
[4] such procedures could take very long to determine whether or not two graphs of order 15 are 
isomorphic. Hence, such methods are clearly inadequate for large order graphs. 

If a computer is to be used to determine whether or not two graphs are isomorphic, the graphs must 
be represented in a form suitable for machine processing. The adjacency matrix of a graph is such a 
standard form. First, the nodes of the graph are arbitrarily associated with the integers 1,2,3 . . . n where 
n is the order of the graph. The (i 9 j) component of the adjacency matrix is then defined to be 1 if there 
is an edge connecting the ith node and the yth node and is otherwise. The (i 9 i) entry commonly con- 
tains the label given the node. In the case of an unlabeled graph, (i 9 i) entries are typically zero. While 
different orderings of the nodes produce different adjacency matrices, two adjacency matrices are 
isomorphic if and only if, independent of labeling, they correspond to the same graph. Figure 2 
illustrates this. 




Figure 2a 



Figure 2b 



A 2 = P^P Where P 



/o 1 o\ 
10 
1 
10 
10 0/ 



l 



Figure 2c 



The algorithm about to be described is divided roughly into three sections each satisfying a specific 
purpose. The first section consists of several basic tests for nonisomorphism. These tests are all 0(n s ) or 
better and are surprisingly effective in distinguishing nonisomorphic graphs. Since nearly all pairs of 
graphs that pass through the initial tests undistinguished are likely to be isomorphic, the second part of 
the algorithm involves the attempted construction of a permutation between the adjacency matrices of 
the graphs. In the great majority of cases, this is easily and quickly accomplished. The second section 
simultaneously serves as a more complicated test for nonisomorphism should it be found impossible to 
generate a satisfactory permutation. 



1 Figures in brackets indicate the literature references at the end of this paper. 
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The final section is specifically designed to deal with certain classes of matrices that do not easily 
yield to the tests or efforts to construct a permutation found in the previous two components. In general 
the tests in this section are more time consuming, generally 0(ra 4 ), and are quite detailed. Methods of 
finding a satisfactory permutation are included as well as tests for nonisomorphism. 

A complete FORTRAN listing of the algorithm, including explanatory comments, follows its 
description in section 3 and the test results provided in section 4 demonstrate the efficiency of the 
algorithm. 

2. The Algorithm 

Initially, the algorithm employs four basic tests, each of which checks for the failure of a condition 
necessary for isomorphism. The first test determines whether or not the node labels of the two graphs 
are identical. If they are not, then the corresponding graphs are not isomorphic. In the case of unlabeled 
graphs, the test yields no information and is unnecessary. This test is performed by SUBROUTINE 
LABEL and takes 0(n 2 ) amount of time when n is the order of the graphs. 

A more useful test involves the comparison of the node valences of the two graphs. The valence of 
a node is the number of edges connected to it. The valence of each node, once computed, is combined 
with the node's label to form a new label for that node. Two nodes of a graph will have the same new 
label if and only if they have the same valence as well as identical initial labels. Thus two graphs with 
different sets of new labels are themselves nonisomorphic. For example, the graphs in figure 3 share the 
same set of node valences and of initial labels yet when these two qualities are considered mutually, the 
graphs are found to be nonisomorphic. 
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Figure 3 



The set of initial node labels is (A,A,B,C,C,C,C) for each graph. Similarly the set of node valences 
is (1,1,1,2,2,3,4) in each case. When combined to produce sets of new labels, however, the following 
two distinct sets are formed: (A3,A4,B2,C1,C1,C1,C2) and (A2,A3,B4,C1,C1,C1,C2) . Since these two 
sets are clearly not identical, the two graphs are not isomorphic. This test takes 0(n 2 ) time and is per- 
formed in SUBROUTINE VALENC. 

The third basic test involves the search for and labeling of duplicate nodes. A given node is called 
an ordinary duplicate of degree k if there exist k — 1 and no more than k — 1 other identically labeled 
nodes in the graph which are connected by edges to the same nodes as is the given node. A given node 
is called a connected duplicate of degree A: if the A: — 1 nodes just mentioned are all mutually connected. 
To further illustrate consider the graph in figure 4. The numbers are only for use in referencing the 
nodes while the letters correspond to actual labels. 

Nodes 1 and 3 are connected duplicates of degree 2. Nodes 2 and 4 are ordinary duplicates of 
degree 2. Nodes 6 and 7 are not duplicates since they have different labels. Finally nodes 10, 11, and 12 
are ordinary duplicates of degree 3. All other nodes are trivially duplicates of degree 1. 

Node duplicate degrees are easily found from the 0,1 adjacency matrix of a graph. Groups of 
ordinary duplicate nodes correspond to groups of identical rows (with identical labels) of the adjacency 
matrix with 0's entered on the diagonal. If l's are entered on the diagonal, then the connected dupli- 
cates can be similarly found. It is easily seen that two graphs are not isomorphic if the ordinary and 
connected duplicate degrees of their nodes differ in any manner. 
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Figure 4 



As before, the node labels of each graph are updated to reflect the newly obtained information and 
then compared. This is done in SUBROUTINE DUPLIC and takes 0(rc 3 ) time. 

While the identification of duplicate nodes serves as a worthwhile test for nonisomorphism, a 
potentially much more useful purpose is also served. A graph with a large number of nodes with 
duplicate degrees greater than one yields an adjacency matrix, A 7 and a large number of nontrivial 
automorphisms, that is permutations P such that P T AP = A. For reasons explained further on, graphs 
whose adjacency matrices have this property are especially difficult to analyze and any procedure which 
efficiently reduces the number of such permutations through appropriate labeling of the nodes of the 
graph, will greatly reduce the complexity of following computations. Labeling of duplicate nodes was 
found to be a particularly fast and simple way to achieve this end for many graphs. It should be noted 
that this labeling alone does not necessarily reduce the automorphism group to the identity. 

A well-known necessary condition for isomorphism involves the comparison of spectra (set of 
eigenvalues) of the 0,1 adjacency matrices of two graphs. Much effort has been devoted to nonisomor- 
phic pairs for which the 0,1 adjacency matrices have the same spectrum. Such pairs of graphs are called 
cospectral [5]. Many cospectral pairs of graphs are distinguishable, however, through comparison of 
the spectra generated from adjacency matrices with numbers different from 0,1 used to represent non- 
edges and edges [6]. It appears that the use of and 1 in the computation of the spectra of graphs 
ignores much of their inherent structure. A greater amount of information about the graph can be 
discovered from the examination of the spectra of adjacency matrices with l's representing no edges, 
and a variable x representing an edge. Due to the nature of the characteristic polynomial, no generality 
is lost by replacing x with either a transcendental number or certain large functions of the number of 
nodes in the graph. Unfortunately neither option is feasible for use with the computer since a trans- 
cendental number is truncated to a rational and a large value of x yields errors in computation small in 
comparison to x but large in comparison to 1. A compromise solution was developed that allows x to be 
a function that increases linearly with the order of the graph yet is essentially as effective in distin- 
guishing nonisomorphic graphs as is the more general procedure. 

In order to utilize the information stored in the updated labels of the nodes, the diagonal elements 
of the adjacency matrices are assigned the value of the appropriate label. The eigenvalues of each 
modified adjacency matrix are then computed and compared. Should the resulting spectra not be 
identical, the graphs are not isomorphic. The preparation of the adjacency matrices for eigenvalue 
computation is effected in SUBROUTINE PREPAR and takes 0(n 2 ) time. The computation of the 
eigenvalues as well as the eigenvectors (for future use) of each matrix is carried out in SUBROUTINE 
SPECTR and takes 0(n 3 ) time. 

Since two graphs that are not distinguished by the preceding sequence of tests may well be isomor- 
phic, it is reasonable at this point to attempt to generate a permutation that relates the two graphs. The 
number of such possible permutations can be greatly limited with efficient use of the information stored 
in the eigenvectors of the two adjacency matrices. If A and B are two symmetric matrices with identical 
spectra then there exist orthogonal matrices U and V and a diagonal matrix D such that 

U T AU = D and V T BV = D. 
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Further assuming that A and B are permutationally equivalent we have B = P T AP where P is a 
permutation matrix. Then : 

JJ T AV— V T P T APV 

D— V T P T UDU T PV 

U T PVD = DU T PV 

and thus D commutes with E = U T PV. 



If E — {eij) and D = (c? i; -) with d. u = for i ^ /*, then Z)jF = (dney) and £/) = {d^eij). Thus if 



^< 7^ ^ tnen c?>y^ii = rf^y — > Cy = 



Thus E : 



£ 2 











where Ei is an orthogonal matrix of degree equal to the multiplicity of the ith distinct eigenvalue of A. 

The relation UE = PV serves to greatly restrict the number of permutations, P. Consider first the 
restrictions imposed by the existence of eigenvalues of multiplicity one. In this case, the corresponding 
matrix E L is trivially (1) or ( — 1). Thus for each / such that djj is unique, we know that the set of 
permutations of 



combined with 



r »w i 
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those of 



all permutations of A into B. 
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V 2 j 
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< Vnj > 



are simply the 



eigenvectors of A and B corresponding to the eigenvalue d n and are easily computed. 

Since the components of the eigenvectors of most modified adjacency matrices are unique, Jie 
number of possible permutations of A into B is small. In these cases, all possible permutations are 
quickly generated and checked until one is found that relates A and B or until the supply is exhausted. 
In the latter case, the graphs are not isomorphic. 

As was mentioned earlier graphs that yield an adjacency matrix, A, and several permutations, P, 
such that P T AP = A are very difficult to work with. The reason for this phenomenon now becomes 
apparent. If the number of permutations of A into itself is large, then there will be many pairs of k and 
/ such that Ukj = Uij for every j such that djj is unique. This, in turn, greatly increases the number of 
permutations 
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of 



U 2 j 



lifij j 



into 



vis 

v 2 j 



Vnj 



and thus 



the amount of time which may be required to find a permutation of A into B and certainly the amount 
of time required to exhaust all possibilities is greatly increased. As an example, certain stochastic 
adjacency matrices have only one eigenvalue of multiplicity one. In many of these cases, the corresponding 
eigenvector components are all identical. In this, as well as in other less severe cases, the above described 
method of finding a permutation of A into B is not practical. 

For such graphs a more sophisticated, as well as time consuming, approach is required. Again 
referring to the relation UE = PV, it is useful to examine the case when an eigenvalue has multiplicity 
two. In other words, assume that djj = d kk for k =/, / + 1 but djj y^ d^k for all other k. Then we know 
that 



U 2 j U 2 j+1 



* U n j U n j+\ J 



U\j l^iy+i 
U 2 j #2/+l 



t, U n j U n j + i 



(;-!)=■ 



or that 



(: -I) 



vy i>i/+i 

V 2 j ^2j+l 



V n j V n j + i 



V\j v lj+1 
V 2 j v 2j+1 



<m Vnj Vnj+1 * 



In each case, there are n possible values for a. Selecting i such that either u^ =^ or u ij+1 ^ and 
assuming that there is a permutation of A into B which takes the ith node of A into the &th node of B, 
we have 



Uijd + u ij+ ib 


— V k j 


Uijb -f- Uij+id 


= Vjcj+l 


or 




Uijd -f u ij+1 b 


= v kj 


Uijb — Uij+id 


= Vlcj+1 



which yields 



/ u^ Ui j+1 \ ( a \ __ / v kj \ 

\ Ui j+1 —u^ J \ b ) \ v kj+1 ) 



or 



/ uij u ij+1 \ ( a \ __ / vy \ 

\ -~Ui j+1 u^ / \ b ) \ v kj+1 J 



Solving for a and b yields : 



a = 



Ujj Vkj "I" Uij+1 v kj+l 

Uij 2 + u ij+1 2 



» — Uij Vkj+l ~f" Uij+1 v kj 

Uij 2 + u ij+1 2 
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or 



b = 



Ujj Vkj — Ujj+i Vkj+1 
2 I 2 

Mty Vkj+1 ~T~ Ujj+1 v kj 



Uij 2 + Ity+i 2 

Note that a 2 + & 2 = 1 so that \a\ < 1. Thus the permutations P are limited to those that satisfy 



U 2 j 



+ b 



Ulj+i 

U>2j+l 



v. U n j + i 



= P 






for any of the In pairs of {a«b) , two pairs corresponding to each value of k. For each (a,6), the prob- 
lem then reduces to the previously examined case of eigenvalues of multiplicity one. The only difference 
is that the process must be repeated 2n times before all possible permutations have been checked. 

SUBROUTINE OPTEIG selects an optimal eigenvector, one that has as few multiple components 
as possible, and takes 0(ti 3 ) time. Tf the selected eigenvector corresponds to an eigenvalue of multi- 
plicity one, SUBROUTINE SINGLE is called upon. This subroutine is 0(/i 3 ) but the leading coefficient 
can be very large if there are large numbers of multiple eigenvector components. SUBROUTINE 
DOUBLE calls SUBROUTINE SINGLE at most In times and thus is 0(n 4 ) with the possibility of a very 
large leading coefficient. 

Due to the complex nature of the relation UE = PV for blocks of E of dimension greater than two, 
little progress was made in obtaining meaningful data from the eigenvectors corresponding to eigen- 
values of multiplicity greater than two. 

Should neither the eigenvectors corresponding to eigenvalues of multiplicity one nor those of multi- 
plicity two sufficiently restrict the number of permutations possible, still other algebraic approaches can 
be employed. 

One such approach is the examination of the spectra of the n subgraphs of order n — 1 of the 
original graph. An obvious necessary condition for isomorphism is that the two graphs generate the 
same sets of spectra. Further there is the possibility that in a given set of spectra, many of them will be 
unique, thus creating the possibility of generating and testing all possible permutations in a small 
amount of time. Distinguishing node duplicates again helps to reduce the maximum number of permu- 
tations possible. Experience with a limited number of large stochastic matrices has shown, however, that 
this procedure is far more effective as a nonisomorphism test than as one that attempts to generate and 
test all possible permutations. This procedure is contained in SUBROUTINE SUBVAL and takes 0(rc 4 ) 
time. 

Carrying the eigenvector analysis further all possible permutations can be derived through exam- 
ination of eigenvectors of subgraphs of the original graph. This method has proved to be one of the 
most successful in distinguishing very difficult large, stochastic graphs. For simplicity and speed, all 
submatrices of order n — 1 are searched for the optimum or a sufficiently satisfactory eigenvector 
corresponding to an eigenvalue of multiplicity one. This eigenvector is then compared with the cor- 
responding eigenvector generated from each of the subgraphs of the second graph that has a spectrum 
matching the spectrum of the subgraph which generated the optimal eigenvector. This procedure is 
executed in SUBROUTINE SUBVEC. The order is 0(/i 4 ) as SUBROUTINE SINGLE can be called a 
maximum of n times. This approach tends to be more time consuming than the others mentioned since 
it could require the generation of n eigenvector matrices. 

In a similar fashion, one could examine the eigenvectors of subgraphs corresponding to eigenvalues 
of multiplicity two. Another possible approach would be to examine the n 2 spectra generated by 
sequentially generating all possible subgraphs of order n — 2. These and other possible procedures are 
0(n 5 ) or worse and very time consuming. It does not appear that the trade off of speed for effectiveness 
dictates inclusion of such unwieldly approaches in the algorithm. 



3. Computer Implementation 

The algorithm was organized into fourteen subroutines and coded in FORTRAN V for use with 
the UNIVAC 1108 computer. SUBROUTINE GISOM serves as the coordinating subroutine and is the 
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only subroutine that need be referenced by a main program. Each subroutine listing is fully commented, 
indicating values of variables on input, output and during the life of the subroutine. 

All eigenvalues and eigenvectors were calculated using the EISPACK software described in [9]. 
Reference to [9] should be made for a detailed explanation of the function, running time, error mes- 
sages and accuracy of the EISPACK subroutines. Throughout the program two numbers, A and B are 
considered identical ii \\A\ — \B\\ < MACHP * (1 -f \A\) . The precision constant, MACHP, is deter- 
mined and set by the user in the call to GISOM. For use with the EISPACK software on the UNI VAC 
1108, a value of MACHP = 0.0003 was found to be satisfactory. 

Another user option in the call to GISOM command sets a limit on the number of permutations 
that may be tested for a given pair of graphs. Generally, this value, MAXT, can be set to a very large 
number, say 5000, since very little time is expended on each test. MAXT is of greatest use when it is 
desired to sequentially test several pairs of very large matrices without spending too much time on 
any one pair. 

Time rather than storage conservation was emphasized. The subroutines utilize five n X n arrays 
and several single dimension 0(n) arrays, all of which are passed as part of the subroutine call state- 
ments. There are no common blocks and no input-output statements in the subroutine as all pertinent 
data is returned through the CALL GISOM command in the main program. The STATUS and ERR 
variables serve as the key to what happened when the two graphs were tested for isomorphism. The 
other variables contain a variety of data whose meaning depends on the value of ERR and STATUS. 
This data can often be quite useful in determining exactly why two graphs are not isomorphic. The 
interpretation of ERR and STATUS is given with the commentary of each subroutine. A list of the 
program follows. 
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SUBROUTINE GISOM C A , B , C ,XA , XB , DA , DB , EA , EB , FA , FB , SC , SD , S£ , CE , CF , G , H 
+ , P , STATU S , ERR , MACHP , MAXT , N , NM , NT ) 

C 

C THIS SUBROUTINE IS THE CONTROL SUBROUTINE FOR THE ALGORITHM. 
C 

C ON INPUT: 

C A AND B CONTAIN THE 0,1 FORM OF THE ADJACENCY MATRICES OF THE TWO 
C GRAPHS TO BE COMPARED. THE DIAGONAL ELEMENTS CONTAIN LABELS 
C WHICH MAY RANGE IN VALUE FROM TO 99. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. A 
C VALUE OF .0003 IS SUGGESTED WHEN THE PROGRAM IS RUN ON A 
C MACHINE COMPARABLE TO THE UNIVAC 1108 AND WITH SOFTWARE 
C EQUIVALENT TO EISPACK IN ACCURACY. 

C MAXT IS APPROXIMATELY THE MAXIMUM NUMBER OF PERMUTATIONS THAT MAY 
C BE TESTED. EXCEPT IN CASES OF SEVERE TIME RESTRICTIONS, A 
C VALUE OF 500 IS SUGGESTED. 

C N IS THE ORDER OF THE GRAPHS TO BE COMPARED. 
C NM IS THE DIMENSION OF MOST OF THE ARRAYS. 
C NT = 2*NM AND IS THE DIMENSION OF CE AND CF . 
C 

C ON RETURN: 

C A AND B ARE POSSIBLY DESTROYED. 

C P CONTAINS A PERMUTATION SUCH THAT AC I , J ) =B ( P ( I ) , P ( U ) ) FOR ALL I. , J 
C SHOULD SUCH A PERMUTATION BE FOUND TO EXIST. 

C STATUS RANGES IN VALUE FROM 1 TO 1 AND RECORDS THE POINT IN GISOM 
C WHERE THE ALGORITHM TERMINATED. FOR EXAMPLE, A VALUE OF STATUS 
C =3 MEANS THAT THE ALGORITHM TERMINATED FOLLOWING THE NODE 
C VALENCE CHECK WHICH IMPLIES THAT THE NODE VALENCES OF THE TWO 
C GRAPHS DIFFERED AND THUS THAT THE TWO GRAPHS ARE NOT 
C ISOMORPHIC. 

C ERR CONTAINS INFORMATION PERTAINING TO THE REASON FOR THE 
C TERMINATION OF THE ALGORITHM. AN EXPLANATION OF ITS MEANING 
C CAN BE FOUND WITH THE COMMENTS IN THE SUBROUTINE THAT WAS LAST 
C CALLED BY GISOM. FOR EXAMPLE, A STATUS VALUE OF 7 COMBINED 
C WITH AN ERR OF IMPLIES THAT A PERMUTATION WAS FOUND THAT 
C RELATES A AND 3 THROUGH USE OF AN EIGENVECTOR CORRESPONDING TO 
C AN EIGENVALUE OF MULTIPLICITY ONE. 
C MACHP, MAXT, N, NM, AND NT ARE UNCHANGED. 
C 

C ALL OTHER VARIABLES ARE CLASSIFIED AS WORKING VARIABLES. SUCH 
C VARIAELES MAY CONTAIN RELEVANT INFORMATION DEPENDING ON THE VALUES OF 
C STATUS AND ERR. 
C 

INTEGER N,NM,NT , ERR , MAXT , G ( NM ) , H ( NM ) , PC NM ) , STATUS, UA 

REAL A(NM,NM) ,BCNM,NM) ,CCNM,NM) ,XACNM,NM) ,XBCNM,NM) ,DACNM) , D3CNM) , 
+ EACNM) ,EBCNM) ,FACNM) , FB( NM) , SC ( NM) , SDC NM ) ,SECNM) ,CECNT) , MACHP 
+ ,CFCNT) 
C 

ERR=0 

STATUS=1 

IF CN.GT.NM) RETURN 
C 

C COMPARE THE NODE LABELS OF THE TWO GRAPHS. 
C 

STATUS =2 

CALL LABEL C A , B , DA , D3 , SC , ERR , N , NM ) 

IF CERR.GT.0) RETURN 
C 
C COMPARE THE NODE VALENCES AND UPDATE THE CORRESPONDING NODE LABELS. 
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c 

STATU S=3 

CALL VALENC C A , 5 , DA , DB , SC , ERR , N , NM) 

IF (ERR.GT.O) RETURN 
C 

C FIND, LABEL, AND COMPARE ORDINARY AND CONNECTED DUPLICATE NODES. 
C 

1 STATUS=STATUS+1 

CALL DUPLIC (A,B,DA,DB, SC, SD, ERR, STATUS, N,NM) 
IF (ERR.GT.O) RETURN 

IF (STATUS. EQ. 4) GO TO 10 
C 

C CONVERT THE 0,1 FORM OF THE ADJACENCY MATRIX TO THE 1,X FORM. 
C 

CALL PREPAR ( A , B , DA , DB , N , NM ) 
C 

C COMPUTE AND COMPARE THE SPECTRA OF THE TWO ADJACENCY MATRICES. 
C 

STATUS=6 

CALL SPECTR ( A , B , XA , XB , EA , EB , SC , G , ERR , MACHP , N , NM ) 

IF (ERR.NE.O) RETURN 
C 

C LOCATE A SUITABLE OR OPTIMAL EIGENVECTOR. 
C 

CALL OPTEIG (XB , DB, G , ERR , MACHP , JA , N , NM) 

GO TO (20 ,30 ,50 ) , ERR 
C 

C ATTEMPT TO FIND A PERMUTATION THAT RELATES A AND 3 THROUGH USE OF AN 
C EIGENVECTOR THAT CORRESPONDS TO AN EIGENVALUE OF MULTIPLICITY ONE. 
C 

20 STATUS=7 

CALL SINGLE ( A , B ,XA ,X3 , FA , FB , H , P , ERR , MACHP , MAXT , JA, N , NM ) 

GO TO 40 
C 

C ATTEMPT TO FIND A PERMUTATION THAT RELATES A AND B THROUGH USE OF AN 
C EIGENVECTOR THAT CORRESPONDS TO AN EIGENVALUE OF MULTIPLICITY TWO. 
C 

30 STATUS=8 

CALL DOUBLE ( A , 3 , C , XA , XB , CE , CF , FA , F3 , H , P , ERR , MAXT , MACHP , JA , N , NM , NT 
+ ) 

40 IF (ERR.NE.4) RETURN 
C 

C COMPUTE AND COMPARE THE SPECTRA OF THE N ORDER N- 1 SUBGRAPHS OF THE 
C ORIGINAL GRAPHS. 
C 

5 STATUS=9 

CALL SUBVAL ( A , B , C ,XA, SC , SD , SE ,H , P , ERR , MACHP , MAXT , N , NM ) 

IF (ERR.NE. .AND.ERR.NE.5+N) RETURN 
C 

C ATTEMPT TO FIND A PERMUTATION THAT RELATES A AND B THROUGH USE OF AN 
C EIGENVECTOR THAT CORRESPONDS TO AN EIGENVALUE OF MULTIPLICITY ONE IN 
C A SUBGRAPH OF ORDER N- 1 . 
C 

STATUS* 1 

CALL SUBVEC ( A , B , C ,XA ,XB , DA , DB , EA , EB, FA , F3 , SC , SD, SE , G , H ,P , ERR , 
+ MACHP, MAXT, N,NM) 

C 

RETURN 

END 
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SUBROUTINE LABEL ( A , 3 , DA , DB , SC , ERR, N , NM) 
C 

C THIS SUBROUTINE RECORDS AND COMPARES THE INITIAL LABELS OF THE NODES 
C OF THE TWO GRAPHS TO BE COMPARED. 
C 

C ON INPUT: 

C A AND B CONTAIN THE TWO ADJACENCY MATRICES TO BE COMPARED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A AND B ARE UNCHANGED. 

C DA AND DB CONTAIN THE DIAGONAL ELEMENTS -LABELS- OF A AND B. 
C ERR = IF DA AND DB ARE ISOMORPHIC. OTHERWISE ERR IS THE INDEX 
C OF THE ELEMENT OF DA WHICH CANNOT BE MATCHED IN DB. 
C N AND NM ARE UNCHANGED. 
C 

C WORKING ARRAYS : SC 
C 

INTEGER N,NM,ERR, I 

REAL ACNM,NM) ,B(NM,NM) , DAC NM) , DBC NM ) , SC ( NM) 
C 

C RECORD THE LABELS. 
C 

DO 10 I = 1,N 

DAC I )=A( 1,1) 
DBC I) = B( 1,1) 
10 CONTINUE 
C 

C COMPARE THE LABELS. 

C NOTE THAT SINCE THE LABELS ARE INTEGERS, THE MACHINE PRECISION 
C CONSTANT IS ZERO. 
C 

CALL COMP C DA,DB,SC,ERR, ,N,NM) 
C 

RETURN 

END 
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SUBROUTINE COMP C HA , RB , RC , ERR ,MACHP , N ,NM > 
C 

C THIS SUBROUTINE DETERMINES WHETHER OR NOT TWO VECTOR ARRAYS ARE 
C ISOMORPHIC* 
C 

C ON INPUT: 

C RA AND RB ARE THE TWO VECTOR ARRAYS TO BE COMPARED. 
C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO 5E TESTED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C RA AND RB ARE UNCHANGED. 

C ERR = IF RA AND RB ARE ISOMORPHIC. OTHERWISE ERR IS THE INDEX 
C OF THE ELEMENT OF RA WHICH CANNOT BE MATCHED IN RB. 
C MACHP, N, AND NM ARE UNCHANGED. 
C 

C WORKING ARRAYS: RC 
C 

INTEGER NM,N,ERR,I,J 

REAL RA(NM) , RBCNM) ,RC(NM) , MACHP, VAL 
C 

VAL=1 
ERR=0 
C 

C COPY RB INTO RC AND COMPUTE A RELATIVELY LARGE NUMBER. 
C 

DO 10 I=1,N 

RCC I)=RB( I) 
VAL=VAL+ABS(R3C I) ) 
10 CONTINUE 
C 

C COMPARE THE TWO ARRAYS. 
C 

DO 3 I = 1 , N 

DO 20 J=l ,N 

IF CABSCRAC I )-RCC J) ) .GT.MACHP*< i+ASSCRCC J) > )> GO TO 20 
RC(U)=VAL 
GO TO 3 
2 CONTINUE 
ERR=I 
RETURN 
30 CONTINUE 



C 



RETURN 
END 
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SUBROUTINE VALENC ( A , B , DA , DE , SC , ERR , N , NM) 
C 

C THIS SUBROUTINE COMPUTES AND COMPARES THE NODE VALENCES OF THE 

C TWO GRAPHS. 
C 

C ON INPUT: 

C A AND B CONTAIN THE TWO ADJACENCY MATRICES TO BE COMPARED. 
C DA AND DB CONTAIN THE NODE LABELS OF A AND B. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON .RETURN: 

C A AND B ARE UNCHANGED. 

C DA AND DB CONTAIN THE UPDATED NODE LABELS. 

C ERR = IF DA AND DB ARE ISOMORPHIC. OTHERWISE ERR IS THE INDEX 
C OF THE ELEMENT OF DA WHICH CANNOT BE MATCHED IN DB. 
C N AND NM ARE UNCHANGED. 
C 

C WORKING ARRAYS: SC 

c ! 

INTEGER NM,N,ERR, I , J i 

REAL ACNM,NM) ,B(NM,NM) , DA<NM) ,DB(NM) , SC(NM) 
C 

DO 20 1 = 1 ,N 
C 

C MAKE ROOM FOR NEW INFORMATION. 
C 

DA (I )=10 0*DACI) 
DEC I ) = 100*DE( I) 
C 

C COMPUTE AND RECORD NODE VALENCES - UPDATE LABELS. 
C 

DO 10 J=l ,N 

IF C I.EQ. J) GO TO 10 
DA( I )=DA(I)+A( I , J) 
DBC I )=DB( I )+E( I ,U) 
10 CONTINUE 
2 CONTINUE 
C 

C COMPARE UPDATED LABELS. 
C 

CALL COMP (DA,DB,SC,ERR,0 ,N,NM) 
C 

RETURN 
END 
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SUBROUTINE DUPLIC < A, E, DA , DB, SC , SD , ERR , Q , N, NM) 
C 

C THIS SUBROUTINE LOCATES AND DISTINGUISHES DUPLICATE NODES THROUGH 

C LABEL MODIFICATION, 
C 

C ON INPUT: 

C A AND B CONTAIN THE TWO ADJACENCY MATRICES TO BE COMPARED. 
C DA AND DB CONTAIN CURRENT NODE LABEL DATA. 

C Q = 4 IF O'S ARE TO BE PLACED ON THE DIAGONAL OF A AND B, AND 
C Q = 5 IF l'S ARE TO BE USED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A AND B HAVE DIAGONAL ENTRIES OF ALL O'S OR ALL l'S BUT ARE 
C OTHERWISE UNCHANGED. 

C DA AND DB CONTAIN UPDATED LABEL DATA. 

C ERR « IF DA AND DB ARE ISOMORPHIC. OTHERWISE ERR IS THE INDEX 
C OF THE ELEMENT OF DA WHICH CANNOT BE MATCHED IN DB. 
C Q, N, AND NM ARE UNCHANGED. 
C 

C WORKING VARIABLES: SC,SD 
C 

INTEGER NM,N,NN, I ,U,K,L,ERR,Q 

REAL A(NM,NM) , BC NM , NM) , DACNM) ,D3(NM) ,SC(NM) ,SD(NM) 



NN=N-1 

DO 10 1 = 1, N 



C 

C SET DIAGONAL ELEMENTS TO OH 1. 

C 

AC I, I) = Q-4 

B(I, I)=Q-4 
C 

C MAKE ROOM FOR NEW DATA. 
C 

DACI)=100*DACI) 

DB(I>=10G*DB(I ) 
C 

C STORE DA IN SC AND DB IN SD. 
C 



SCCI)=DA(I) 

SD(I)=DBC I ) 
10 CONTINUE 
C 

DO 6 1=1 ,NN 
L=I + 1 

DO 5 K=L,N 
C 

C SEARCH FOR IDENTICAL ROWS IN A. 
C 

IF CSCC I).NE.SC(K)) GO TO 30 
DO 20 J=1,N 

IF CACI ,J) .NE.A(K,J) ) GO TO 30 
20 CONTINUE 

DA(K)=DA(K)+1 
DAC I)=DA( I )+2 
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c 

C SEARCH FOR IDENTICAL ROWS IN B. 
C 

30 IF (SDC I ).NE. SDCK)) GO TO 50 

DO 40 J«l ,N 

IF (EC I ,J) .NE.BCK, J) ) GO TO 53 
40 CONTINUE 

DB(K)=DB(K)+1 
DB< I )=DB( I )+2 
5 CONTINUE 
60 CONTINUE 
C 

C COMPARE UPDATED LABELS. 
C 

CALL COMP (DA,DB,SC,ERR,Q ,N,NM) 
C 

RETURN 
END 
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SUBROUTINE PREPAR ( A , B , DA , DB , N , NM) 
C 

C THIS SUBROUTINE CONVERTS THE ADJACENCY MATRIX FROM 0,1 FORMAT TO 1,X 
C FORMAT. IN ADDITION THE LABELS ARE CONVERTED INTO MULTIPLES OF N 
C AND ARE ENTERED AS DIAGONAL ELEMENTS OF A AND B. 
C 

C ON INPUT: 

C A AND 3 ARE THE ADJACENCY MATRICES IN , 1 FORMAT, 
C DA AND DB CONTAIN CURRENT NODE LABEL DATA. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A AND B CONTAIN* THE TWO ADJACENCY MATRICES IN 1,X FORM WITH 
C MODIFIED DIAGONAL ELEMENTS IN MULTIPLES OF N. 

C DA AND DB CONTAIN THE NODE LABEL DATA IN THE FORM OF MULTIPLES 
C OF N. 

C N AND NM ARE UNCHANGED. 
C 

INTEGER NM,N,L 

REAL A(NM,NM) ,B(NM,MM) ,DACNM) ,DB(NM) , FN , NA , NB , QA , QB, INF 
C 

C FN IS THE FUNCTION FOR X-l IN THE 1,X FORM OF THE ADJACENCY MATRIX. 
C 

FN=N/2 
QA=-1 
Q3=-l 
NA=0 
N3=0 

I NF=1 000000000 
C 

C CONVERT A AND B TO 1,X FORM. 
C 

DO 20 1 = 1, N 

DO 10 J=l ,N 

AC I ,J)=A( I ,J)*FN+1 
BCI ,J)=B( I ,J)*FM+1 
10 CONTINUE 
20 CONTINUE 
C 

DO 5 I = 1 , N 
C 

C LABEL DIAGONAL ELEMENTS OF A IN MULTIPLES OF W. 
C 

L=l 

DO 3 J=1,N 

IF CDA(J) .LT.DA(D) L = J 
30 CONTINUE 

IF (DA(L).GT.QA) NA=NA+N 
A(L,L)=NA 
QA=DA(L) 
DACL)=INF 
C 

C LABEL DIAGONAL ELEMENTS OF B IN MULTIPLES OF N. 
C 

L=l 

DO 4 J=1,N 

IF (DBCJ).LT.DB(L)) L=J 
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40 CONTINUE 

IF CDBCD.GT.QB) NB=NB+N 

B(L,L)=NB 

QB=DBCL) 

DB(L)=INF 

5 CONTINUE 
C 

C RESTORE DA AND DE. 

DO 6 I = 1 , N 

DA(I)=A(I,I) 
DBC I )=BCI , I) 

6 CONTINUE 

C i 

RETURN 
END 
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SUBROUTINE SPECTR C A, B , XA , XB , EA , EB , SC , G , ERR, MACHP , N , NM ) 
C 

C THIS SUBROUTINE CALCULATES AND COMPARES THE SPECTRA OF THE TWO 
C MODIFIED ADJACENCY MATRICES. IN ADDITION, DATA PERTAINING TO THE 
C MULTIPLICITY OF THE EIGENVALUES IS RECORDED. 
C 

C ON INPUT: 

C A AND B CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS. 
C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A AND B ARE UNCHANGED. 

C XA AND XB CONTAIN THE EIGENVECTORS OF A AND B IN THE SAME ORDER 
C AS THE CORRESPONDING EIGENVALUES. 

C EA AND EB CONTAIN THE EIGENVALUES OF A AND 3 IN NONDECREAS ING 
C ORDER. 

C G OF I IS 1 IF THE I TH AND 1+1 TH EIGENVALUES OF A ARE IDENTICAL 
C AND IS OTHERWISE. 

C ERR IS NEGATIVE IF AN ERROR CONDITION WAS RAISED IN THE EISPACK 
C SOFTWARE. ERR IS IF THE SPECTRA ARE IDENTICAL. OTHERWISE, 
C ERR IS THE INDEX OF THE EIGENVALUE OF A THAT COULD NOT BE 
C MATCHED WITH AN EIGENVALUE OF B. 
C MACHP, N, AND NM ARE UNCHANGED. 
C 

C WORKING VARIABLES: SC 
C 

INTEGER NM, N, ERR, I ,NN, GCNM) 

REAL A(NM,NM) , B(NM,NM) ,XA(NM,NM) ,XB(NM,NM> , EAC NM) , EEC NM ) , SC ( NM) , 
+ MACHP 
C 

NN=N-1 
C 

C FIND THE EIGENVALUES AND EIGENVECTORS OF A. 
C 

CALL TRED2 (NM , N , A , EA , SC , XA ) 

CALL TQL2 ( NM , N , EA , SC ,XA , ERR ) 

IF (ERR.NE.O) GO TO 10 
C 

C FIND THE EIGENVALUES AND EIGENVECTORS OF B. 
C 

CALL TRED2 ( NM ,N , B , EB , SC ,XB) 
CALL TQL2 C NM , N ,EB , SC ,X3 , ERR) 
10 ERR=-ERR 

IF (ERR.EQ.G) CALL COMP ( EA , EB, SC , ERR , MACHP , N , NM) 
C 

C RECORD EIGENVALUE MULTIPLICITY DATA. 
C 

IF (ERR.NE.O) RETURN 

DO 20 1=1 ,NN 
G( I)=0 

IF (ABS(EA(I )-EA( 1 + 1 ) ) .LT. MACHP* C 1+ABS(EA( I ) ) ) ) G( I)=l 
2 CONTINUE 

G(N)=0 
C 

RETURN 

END 
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SUBROUTINE OPTEIG (X3,D3,G p ERR, MACHP ,JA,N,NM) 
C 

C THIS SUBROUTINE SELECTS AN OPTIMAL EIGENVECTOR FOR USE IN GENERATING 
C A SERIES OF PERMUTATIONS. 
C 

C ON INPUT: 

C XB IS THE EIGENVECTOR MATRIX OF B. 
C DB CONTAINS THE DIAGONAL ELEMENTS OF B. 
C G CONTAINS EIGENVALUE MULTIPLICITY INFORMATION. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C XB IS ESSENTIALLY UNCHANGED. THE POSSIBILITY EXISTS THAT ANY TWO 
C EIGENVECTORS CORRESPONDING TO THE SAME EIGENVALUE COULD BE 
C INTERCHANGED. 

C DB AND G ARE UNCHANGED. 

C ERR a 1 IF THE EIGENVECTOR SELECTED CORRESPONDS TO AN EIGENVALUE 
C OF MULTIPLICITY 1. ERR«2 IF THE CORRESPONDING EIGENVALUE HAS 
C MULTIPLICITY 2. £RR=3 IF NO SUITABLE EIGENVECTOR COULD BE 
C FOUND. 
C MACHP IS UNCHANGED. 

C JA IS THE COLUMN INDEX OF THE OPTIMAL EIGENVECTOR, IF ONE EXISTS. 
C N AND NM ARE UNCHANGED. 
C 

INTEGER NM,N, I , J,PN,LT , PX , ERR , UA , L I , LM , L , NN , G ( NM) 
REAL XB(MM,NM) ,DB(NM) , MACHP 
C 

NN=N-1 
ERR = 3 
PN=1 00 
C 

C ALL EIGENVECTORS WITH THREE OR MORE IDENTICAL COMPONENTS OR THAT 
C CORRESPOND TO EIGENVALUES OF MULTIPLICITY THREE OR MORE ARE 
C ELIMINATED FROM CONSIDERATION. EIGENVECTORS THAT CORRESPOND TO 
C EIGENVALUES OF MULTIPLICITY ONE ARE CONSIDERED PREFERENTIALLY TO 
C THOSE WHICH HAVE MULTIPLICITY TWO. EIGENVECTORS WHICH HAVE FEW PAIRS 
C OF IDENTICAL COMPONENTS ARE FAVORED OVER THOSE WITH MANY SUCH PAIRS. 
C 

DO 8 J = 1 , N 
LT = 
PX=1 
C 

C CHECK FOR MULTIPLICITY OF ASSOCIATED EIGENVALUE. 
C 

IF C J. EQ. 1 ) GO TO 20 

IF CGCJ-1 ) .EQ. 0) GO TO 20 
IF (J.EQ.2) GO TO 10 

IF CGCJ-2) .EG. 1 . ) GO TO 80 
10 IF (G(U).EQ. 1 > GO TO 8 

LT=4 
PX = 4 

GO TO 30 
20 IF (GCJ).EQ.O) GO TO 30 

IF CGCU+1 ) .EQ. 1 ) GO TO 80 
LT=4 
PX = 2 
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c 

C CHECK FOR IDENTICAL COMPONENTS. 
C 

3 DO 5 I = 1 , MM 

LI=I+i 

LM=0 

DO 4 L=LI ,N 

IF (ABSCXBC I , J)-XB<L,J)+DB( I )-DBCL) ) . GT.MACH?) GO TO 40 
LT=LT+i 

IF (LM.EQ. 1 ) GO TO 80 
LM=1 

4 CONTINUE 

5 CONTINUE 

IF CLT.GE.PN) GO TO 60 
JA = J 
ERH=PX 

IF (ERH.LT.4) GO TO 7 
ERR = 2 
JA=JA-1 
C 

C INTERCHANGE TWO EIGENVECTORS THAT CORRESPOND TO THE SAME EIGENVALUE 
C SO THAT IN THE RESULTING PAIR THE FIRST WILL BE MORE SUITABLE THAN 
C THE LAST. 
C 

DO 6 I = 1 , N 

TEMP-XBC I , JA) 
XSC I ,UA)=Xe< I , JA+l ) 
XBC I , JA+l )=TSMP 
60 CONTINUE 

70 PN=LT 

C 

C IF EIGENVECTOR IS SUITABLE -NOT NECESSARILY OPTIMAL- RETURN. 
C 

IF (LT. LT.4) RETURN 
80 CONTINUE 
C 

RETURN 
END 
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SUBROUTINE SINGLE ( A , B , XA , XB , FA , FB , H ,P, ERR , MACHP , MAXT , JA , N , NM ) 
C 

C THIS SUBROUTINE GENERATES AND TESTS UP TO MAXT PERMUTATIONS FROM DATA 
C FOUND IN AN EIGENVECTOR CORRESPONDING TO AN EIGENVALUE OF 
C MULTIPLICITY ONE. 
C 

C ON INPUT: 

C A AND B CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS. 
C XA AND XB ARE THE EIGENVECTOR MATRICES OF A AND B THOUGH ONLY THE 
C JA TH COLUMNS OF XA AND XB ARE USED IN THIS SUBROUTINE. 
C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C MAXT IS THE MAXIMUM NUMBER OF PERMUTATIONS THAT CAN BE TESTED. 
C JA IS THE COLUMN OF XA AND XB WHICH CONTAINS THE EIGENVECTORS TO 
C BE EXAMINED. 

C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A, B, XA, AND XB ARE UNCHANGED. 

C P CONTAINS THE LAST PERMUTATION THAT WAS CHECKED AGAINST A AND B. 
C IF ERR=0, THIS PERMUTATION RELATES A AND B. 

C ERR = -1 IF NO SUITABLE PERMUTATIONS COULD BE GENERATED AND THE 
C GRAPHS ARE NOT ISOMORPHIC. ERR= 1 IF THREE IDENTICAL COMPONENTS 
C WERE FOUND IN THE EIGENVECTOR. £RR=3 IF ALL POSSIBLE 
C PERMUTATIONS FAILED WHEN TESTED AGAINST A AND B AND THE GRAPHS 
C ARE NOT ISOMORPHIC. EFR = 4 IF THE MAXIMUM NUMBER OF 
C PERMUTATIONS WERE TRIED WITHOUT SUCCESS. ERR=0 IF A 
C PERMUTATION WAS FOUND TO RELATE A AND B. 
C MACHP, MAXT, JA , N, AND NM ARE UNCHANGED. 
C 

C WORKING VARIABLES: FA,FB,H 
C 

INTEGER NM,N,ERR, JA , I ,M,H(NM) ,P(NM) ,LT,MAXT 

REAL A(NM,NM) ,B(NM,NM) ,XA(NM,NM) ,XB(NM,NM) f FACNM) ,FB(NM) , MACHP 
C 

M = 
C 

C USE THE NODE LABELS TO AIDE IN THE SEPARATION OF IDENTICAL 
C COMPONENTS. 
C 

DO i 1 = 1 , N 

FA( I )=A( I , I )+XA( I , JA) 
FBC I ) = S( I , D+XBCI , JA) 
10 CONTINUE 
C 

C GENERATE A POTENTIAL PERMUTATION AS WELL AS THE DATA NECESSARY TO 
C FIND ALL OTHER POSSIBLE PERMUTATIONS. 
C 

20 CALL PERM ( FA , FB , H , P , ERR , MACHP , LT , N , NM ) 

IF CERR.EQ. 1 ) RETURN 

IF (ERR. EC. -1) GO TO 30 
C 

C TEST UP TO MAXT PERMUTATIONS AGAINST A AND B. 
C 



CALL TEST ( A , B , H , P , ERR ,MAXT , LT , N , NM) 
IF (ERR.NE.3) RETURN 
3 IF CM.EQ. 1 ) RETURN 
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C SHOULD ALL PERMUTATIONS FAIL, REPEAT THE PROCESS WITH ONE EIGENVECTOR 

C NEGATED TO INCLUDE ALL POSSIBLE CASES. 

C 

DO 4 I = 1 , N 

FBC I )=3C I , I)-XB< I , JA) 
4 CONTINUE 
M=l 

GO TO 20 
C 

END 
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SUBROUT I NE PERM ( F A , FB , H , P , ERR , MACHP , LT , N , NM ) 
C 

C THIS SUBROUTINE GENERATES A SINGLE PERMUTATION BETWEEN THE COMPONENTS 
C OF TWO VECTORS AND RECORDS THE DATA NECESSARY TO GENERATE ALL OTHER 
C POSSIBLE PERMUTATIONS BETWEEN THE COMPONENTS OF THE TWO VECTORS. THE 
C VECTORS MUST NOT CONTAIN MOKE THAN TWO IDENTICAL COMPONENTS. 
C 

C ON INPUT: 

C FA AND FB ARE THE TWO VECTORS FROM WHICH ALL PERMUTATIONS WILL 
C BE GENERATED. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 
C FA IS UNCHANGED. 
C FB IS DESTROYED. 

C H CONTAINS THE INDICES OF THE PAIRS OF IDENTICAL COMPONENTS OF FA. 
C P CONTAINS A PERMUTATION OF THE COMPONENTS OF FA INTO THOSE OF FB. 
C ERR = -1 IF THE TWO VECTORS ARE FOUND TO BE NON- I SOMORPH I C. 
C ERR=1 IF THREE OR MORE COMPONENTS OF FA ARE IDENTICAL. 

C ERR = OTHERWISE. 

C MACHP IS UNCHANGED. 

C LT IS THE NUMBER OF IDENTICAL PAIRS OF COMPONENTS FOUND IN FA. 
C N AND NM ARE UNCHANGED. 
C 

INTEGER NM,N , I , J , L , LT ,HCNM) ,PCNM) ,ERR 

REAL FACNM) ,FBCNM), MACHP, INF 
C 

ERR=1 

LT = 

INF=-1 
C 

DO 7 I = 1 , N 
L = -l 
C 

C FIND ALL COMPONENTS OF FB THAT MATCH THE CURRENT COMPONENT OF FA. 
C 



DO 30 J=l ,N 

IF (ABSCFAC I)-FB( J) ).GT. MACHP) GO TO 30 
IF (L) 10,20,80 
10 FB(J)=INF 

PC I )=J 
20 L=L+1 

3 CONTINUE 
C 

IF (L) 40,70,50 

4 ERRa-1 

RETURN 
C 

C RECORD THE INDEX OF THE FIRST HALF OF A IDENTICAL PAIR IN FA. 
C 

50 LT=LT+1 

H<2*LT-1 ) = I 
C 

C FIND AND RECORD THE INDEX OF THE SECOND HALF OF THE IDENTICAL PAIR 
C IN FA. 
C 
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DO 6 J=l ,N 




IF CABSCFAU) 




H<2*LT)=J 




GO TO 7 


60 


CONTINUE 


70 


CONTINUE 




£RR=0 


80 


RETURN 



•FAC J>> .GT.MACHP.OR. I »EQ. J) GO TO 60 



END 
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SUBROUTINE TEST ( A , S , H , P , ERR ,MAXT , LT , N , NM > 
C 

C THIS SUBROUTINE TESTS UP TO MAXT PERMUTATIONS AGAINST THE ORIGINAL 
C GRAPHS. 
C 

C ON INPUT: 

C A AND B CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS. 
C H CONTAINS THE PAIRS OF INDICES OF P WHICH MUST BE INTERCHANGED 
C BEFORE ALL POSSIBLE PERMUTATIONS CAN BEEN TESTED. 
C P CONTAINS THE FIRST PERMUTATION TO BE CHECKED* 

C MAXT IS THE MAXIMUM NUMBER OF PERMUTATIONS WHICH MAY BE TESTED. 
C LT IS THE NUMBER OF PAIRS OF INDICES OF P THAT MAY HAVE TO BE 
C INTERCHANGED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A, B, AND H ARE UNCHANGED. 

C P CONTAINS THE LAST PERMUTATION THAT WAS CHECKED AGAINST A AND B. 
C IF ERR=0, THEN THIS PERMUTATION RELATES A AND B. 

C ERR ■ IF A PERMUTATION WAS SUCCESSFUL IN RELATING A AND B. 
C ERR=3 IF ALL PERMUTATIONS FAILED WHEN TESTED AGAINST A AND B. 
C ERR=4 IF MAXT PERMUTATIONS WERE TRIED WITHOUT SUCCESS. 
C LT, N, AND NM ARE UNCHANGED. 
C 

INTEGER NM,N,LT,ERR,LB, I , J , B IN , MAXT , T , L , K, H ( NM) ,P(NM) ,M 

REAL A(NM,NM) ,BCNM,NM) 
C 

C LB IS THE MAXIMUM NUMBER LESS 1 OF PERMUTATIONS WHICH MAY HAVE TO BE 
C TESTED. 
C 

LB=2**LT- 1 
C 

C BIN CONTROLS WHICH PAIRS OF INDICES ARE TO BE INTERCHANGED EACH TIME. 
C 

SIN = 

T = G 

ERR = 
C 

10 T=T+1 
C 

C TEST THE PERMUTATION IN P AGAINST A AND B. 
C 

DO 3 I = 1 , N 
L=P<I) 
DO 20 U=1,N 
K=PCJ) 

IF (B(L,K) .NE.AC I , J) ) GO TO 4 
20 CONTINUE 
3 CONTINUE 

RETURN 
C 

C CHECK IF ALL POSSIBLE PERMUTATIONS HAVE BEEN TRIED. 
C 

40 IF (BIN.LT.LB) GO TO 50 
ERR=3 
RETURN 
C 
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C CHECK IF MAXT PERMUTATIONS HAVE BEEN TRIED. 
C 

5 IF CT.LT.MAXT) GO TO 6 

ERR = 4 

RETURN 
C 

C GENERATE THE NEXT PERMUTATION TO EE TESTED. 
C 

60 BIN=BIN+1 

c 

C SEQUENTIALLY INTERCHANGE 1 + 1 PAIRS OF VALUES OF P WHERE 2**1 IS THE 
C LARGEST POWER OF 2 THAT DIVIDES BIN. ON THE AVERAGE, TWO PAIRS WILL 
C BE INTERCHANGED EACH TIME, INDEPENDENT OF THE VALUE OF LT . 
C 

M=l 

DO 7 1=1 ,LT 
L=HC2*I-1 ) 
K=HC2*I ) 
J=P(L) 
P(L)=?(K) 
P(K)=J 
M=2*M 

IF CMODCEIN,M) .GT. ) GO TO 10 
70 CONTINUE 
C 

END 
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SUBROUTINE DOUBLE ( A , B , C , XA , X5 , CE , CF , FA , FB , H , P , ERF: , MAXT , MACHP , JF , 

+ N,NM,NT) 

C 

C THIS SUBROUTINE GENERATES AMD TESTS PERMUTATIONS FROM DATA FOUND IN 

C EIGENVECTORS CORRESPONDING TO EIGENVALUES OF MULTIPLICITY TWO. 
C 

C ON INPUT: 

C A AND B CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS. 
C XA AND XB ARE THE EIGENVECTOR MATRICES OF A AND B. 

C MAXT IS THE MAXIMUM NUMBER OF PERMUTATIONS THAT COULD BE TRIED IF 
C THE OPTIMAL EIGENVECTOR CORRESPONDED TO AN EIGENVALUE OF 
C MULTIPLICTY ONE. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C JF AND JF-fl ARE THE TWO COLUMNS OF XA IN WHICH THE PERTINENT 
C EIGENVECTORS ARE STORED. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF MOST OF THE ARRAYS. 
C NT a 2*NM AND IS THE DIMENSION OF CE AND CF. 
C 

C ON RETURN: 

C A, B, XA, AND XB ARE UNCHANGED. 

C P CONTAINS THE LAST PERMUTATION THAT WAS CHECKED AGAINST A AND E. 
C IF ERRaO, THIS PERMUTATION RELATES A AND B. 

C ERR = 4 IF NO PERMUTATION WAS FOUND TO RELATE A AND B BUT NOT ALL. 
C POSSIBLE PERMUTATIONS WERE TESTED DUE TO THE LIMIT IMPOSED 
C BY MAXT AND MAXTN. ERRaO IF A PERMUTATION SUCCESSFULLY RELATED 
C A AND B. ERRa-1,1, OR 3 IF ALL POSSIBLE PERMUTATIONS FAILED 
C WHEN TESTED AGAINST A AND B AND THE GRAPHS ARE NOT ISOMORPHIC. 
C MAXT, MACHP, JA , N, NM , AND NT ARE UNCHANGED. 
C 

C WORKING VARIABLES: C , CE , CF , FA , FB , H 
C 

INTEGER NM,N,ERR,K, I , J , M , JF , NT , L , MAXT , MAXTN , H ( NM) ,P(NM) 

REAL A(NM,NM) ,B(NM,NM) , C(NM,NM) ,XA(NM,NM> ,XB(NM,NM) ,FA(MM) ,FE(NM) , 

+ CE(NT) ,CF(NT) , MACHP , MACHPN , TM 
C 

C RELAX THE REQUIREMENT FOR EQUALITY AS MULTIPLE OPERATIONS ARE TO BE 

C PERFORMED ON THE EIGENVECTOR COMPONENTS. 
C 

MACHPN=N*MACHP/5 
C 

C SET THE MAXIMUM NUMBER OF TESTS THAT MAY BE PERFORMED ON EACH 

C PROPOSED LINEAR COMBINATION OF THE TWO EIGENVECTORS. 
C 



MAXTNa4*MAXT/N 



DO 10 J=1,N 
L = J 

IF (ABSCXAC J, JF) ) . GT . MACHPN.OR. ABSCXAC J, JF+ 1 ) ) . GT. MACHPN) 
+ GO TO 2 

10 CONTINUE 
C 

C COMPUTE ALL POSSIBLE VALUES OF THF. TWO COEFFICIENTS AND RECORD THEM 
C IN CE AND CF. 
C 

2 TM=XA ( L , JF ) **2+XA ( L , JF + 1 ) **2 
DO 5 J=l ,N 
K=2*J 
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CECK-1 )=<XACL, JF)*XB( J, JF)+XACL, JF+1 )*XB< J, JF+1 ) )/TM 
CFCK-1 >=<XACL, JF+1 ) *X8< J , JF ) -XAC L , JF ) *X3 C J , JF + 1 ) )/TM 
CE(K)=CXA(L, JF)*XB< J,JF)-XA(L, JF+i >*XB( J, JF+1 ) )/TM 
CFCK)=CXACL, JF+1 > *XB< J , JF ) +XA< L, JF) *XB< J, JF + i ) )/TM 

C 

C AVOID UNNECESSARY WORK THROUGH DELETION OF DUPLICATED COEFFICIENTS. 
C 

IF C J.EQ. 1 ) CO TO 40 
M = K-2 
DO 3 I = 1 , M 

IF CABSCCEC I) -CECK-1 )) .LT.MACHP) CECK-1 )=2 
IF CABSCCEC I)-CE(K) ) .LT.MACHP) CECK>=2 

3 CONTINUE 

4 IF CABSCCECK) -CECK-1 ) ) . LT.MACHP) CE(K)=2 
C 

5 CONTINUE 
C 

C TEST ALL POSSIBLE COEFFICIENTS. 
C 

L = 

DO 8 J=l ,NT 
C 

C ELIMINATE FROM CONSIDERATION COEFFICIENTS WITH ABSOLUTE VALUE > 1. 
C 

IF CABSCCECJ) ).GT. 1+MACHPN) GO TO 80 
IF CABSCCEC J) ) .GT. 1 ) CECJ)=1 
C 

C PERFORM THE LINEAR COMBINATION TO CREATE A PSEUDO EIGENVECTOR. 
C 

DO 7 I = 1 , N 

CC I ,JF)=CE<J)*XA< I , JF)+CF(J)*XA< I, JF+1) 
7 CONTINUE 
C 

C TREAT THE PSEUDO EIGENVECTOR AS ONE THAT CORRESPONDS TO AN EIGENVALUE 
C OF MULTIPLICITY ONE AND CHECK ALL POSSIBLE PERMUTATIONS. 
C 

CALL SINGLE C A , B , C , X3 , FA , F B , H , P , ERR , MACHPN , MAXTN , JF , N , MM ) 
IF C ERR . EQ . ) RETURN 
IF CERR.EQ.4) L=4 
60 CONTINUE 

IF CL.EQ.4) ERR=4 
C 

RETURN 
END 
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SUBRQUT INE SUBVAL ( A , B , C , XA , SC , SD , S£ , H , ? , ERR, MACHP , MAXT , N , NM) 
C 

C THIS SUBROUTINE GENERATES AND COMPARES THE SPECTRA OF THE ORDER N-l 

C SUBGRAPHS OF THE TWO ORIGINAL GRAPHS. AN ATTEMPT IS MADE TO GENERATE 

C A PERMUTATION BETWEEN THE TWO GRAPHS BASED ON THIS DATA. 

C 

C ON INPUT: 

C A AND E CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 

C MAXT IS THE MAXIMUM NUMBER OF PERMUTATIONS THAT COULD BE TRIED IF 

C THE OPTIMAL EIGENVECTOR CORRESPONDED TO AN EIGENVALUE OF 

C MULTIPLICITY ONE. 

C N IS THE ORDER OF THE GRAPHS. 

C NM IS THE DIMENSION OF THE ARRAYS. 

C 

C ON RETURN: 

C A AND B ARE UNCHANGED. 

C C CONTAINS SPECTRA DUPLICATION INFORMATION. THE I TH ROW OF C 

C CONTAINS THE INDECIES OF ALL NODES IN A WHICH WHEN UNIQUELY 

C MARKED PRODUCE SPECTRA IDENTICAL TO THAT PRODUCED WHEN THE I TH 

C NODE OF A WAS UNIQUELY MARKED. 

C XA CONTAINS THE N SPECTRA GENERATED FROM THE SUBGRAPHS OF A. EACH 

C ROW CONTAINS ONE SPECTRA. 

C SC OF I IS THE NUMBER OF NODES IN A, INCLUDING I, THAT PRODUCED 

C SPECTRA IDENTICAL TO THAT OF I. 

C P CONTAINS A PERMUTATION OF THE SET OF SPECTRA OF SUBGRAPHS OF A 

C INTO THOSE OF B. IF ERR=N+1, THIS IS A PERMUTATION OF A 

C INTO B. 

C ERR IS NEGATIVE IF AN ERROR OCCURS IN THE EISPACK SOFTWARE. IF 

C THE SETS OF SPECTRA ARE NOT ISOMORPHIC, ERR CONTAINS THE INDEX 

C OF THE NODE IN B WHICH WHEN UNIQUELY MARKED, PRODUCED A SPECTRA 

C THAT COULD NOT BE MATCHED WITH A SPECTRA FROM THE SUBGRAPHS OF 

C A. £RR=G IF THE SETS OF SPECTRA WERE ISOMORPHIC BUT THREE OR 

C MORE NODES OF A, WHICH WHEN UNIQUELY MARKED, PRODUCED IDENTICAL 

C SPECTRA. ERR=N+1 IF A PERMUTATION RELATING A AND B WAS FOUND. 

C ERR=N+4 IF ALL POSSIBLE PERMUTATIONS FAILED TO RELATE A AND B. 

C ERR=N+5 IF MAXT PERMUTATIONS WERE GENERATED AND TESTED WITHOUT 

C SUCCESS. 

C MACHP, MAXT, N, AND NM ARE UNCHANGED. 

C 

C WORKING VARIABLES: SD,S£,H 

C 

I NTEGER NM , N , I , J , K , L , ERR , M , MAXT , LT , H ( NM ) , P C MM ) 

REAL ACNM,NM) ,BCNM,NM) ,C(NM,NM) ,XACNM,NM) , SC(NM) , SDCNM) , SECNM) ,FN, 
+ MACHP 
C 

C FN IS THE FUNCTION FOR X- 1 IN THE 1,X FORM OF THE ADJACENCY MATRIX. 
C 

FN=N/2 
C 

C GENERATE THE SPECTRA OF ALL ORDER N-l SUBGRAPHS OF A. 
C 

DO 5 1 = 1 ,N 
PC I ) = 

c 

C COPY A INTO C 
C 



DO 20 K=1,N 
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DO 10 J=l ,N 




C(K,J)=A(K,J) 




10 CONTINUE 




20 CONTINUE 


c 




c 

r 


UNIQUELY MARK THE I TH NODE 


U 


DO 30 J=1,N 




C(J,J)=C(J,J)+C(I ,J) 




C( I ,J)=3*FN 




CC J, I >=3*FN 




3 CONTINUE 




C(I ,I)=CCI,I)+FN 


c 




c 


COMPUTE THE SPECTRA - C IS 


c 





C=A AND ALL NODES CONNECTED TO IT, 



DESTROYED IN THIS COMPUTATION. 

CALL TREDI < NM,N , C , SC , SD, SE ) 

CALL TQL1 CN,SC,SD,ERH) 

ERR=-ERR 

IF (ERR.NE.0) RETURN 

DO 40 J=1,N 

XA( I , J)«SC< J) 

4 CONTINUE 

5 CONTINUE 
C 

C SEQUENTIALLY COMPUTE THE SPECTRA OF THE N- 1 ORDER SUBGRAPHS OF B, 
C CHECKING THAT EACH CAN RE MATCHED WITH A SIMILAR SPECTRA FROM THE 
C SUBGRAPHS OF A. 
C 

DO 11 I = 1 , N 
DO 7 K=1,N 

DO 6 J=l ,N 

C(K,J)=B(K,J) 
60 CONTINUE 

70 CONTINUE 

DO 8 J=l ,N 

C(J,J)=CC J,J)+CCI ,J) 
CCI ,J)=3*FN 
CCJ,I)=3*FN 
80 CONTINUE 

C( I ,1 )=CCI »I )+FN 

CALL TREDI ( NM , N , C , SC , SD , SE ) 

CALL TQL1 (N,SC,SD,ERR> 

ERR=-ERR 

IF (ERR.NE.0) RETURN 

DO 100 J=1,N 

IF (PCJ).GT.O) GO TO 100 
DO 90 K-l ,N 

IF CABSCXAC J,K)-SC(K> ).GT.MACHP*C 1+ABSCSCCK) ) )) 
+ GO TO 100 

90 CONTINUE 

PCU)=I 
GO TO 110 
100 CONTINUE 
ERR=I 
RETURN 
110 CONTINUE 
C 
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LT = 
L=0 
C 

C COMPUTE AND RECORD THE DATA PERTAINING TO SPECTRA DUPLICATION IN A. 
C 

DO 140 I=1,N 
M=0 
DO 130 J=l ,N 

DO 120 K=1,N 

IF (ABSCXAC I , K)-XAC J,K) ) . GT.MACHP*< 1 +ABSCXAC I ,K) ) ) ) 
+ GO TO 130 

120 CONTINUE 

M = M+1 

CC I ,M)=PCJ) 
IF (M.EQ.3) L=l 

IF (L.EQ. l.OR. I.GE. J) GO TO 130 
C 

C RECORD DUPLICATED PAIRS OF SPECTRA WITH THE HOPE THAT THERE ARE NO 
C TRIPLICATES AND SUBROUTINE TEST CAN BE USED TO FIND A PERMUTATION OF 
C A INTO B OR PROVE THAT ONE DOES NOT EXIST. 
C 

LT=LT+1 
H(2*LT-i ) = I 
H<2*LT)=J 
130 CONTINUE 

SCCI)=M 
140 CONTINUE 

IF CL.EQ. 1) RETURN 

CALL TEST ( A , B , H , P , ERR , MAXT , LT , N , NM ) 

ERR=ERR+N+1 



C 



C 



RETURN 
END 
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SUBROUTINE SUBVEC ( A , B , C , XA ,XB, DA , DB , EA , EB, FA, FB, SC , SD, SE , G , H , P , 
+ ERR , MACHP , MAXT , N , NM ) 

C 

C THIS SUBROUTINE ATTEMPTS TO GENERATE A PERMUTATION BETWEEN A AND B 
C THROUGH EXAMINATION OF EIGENVECTORS CORRESPONDING TO EIGENVALUES OF 
C MULTIPLICITY ONE OF ORDER N-l SUBGRAPHS OF A AND B. 
C 

C ON INPUT: 

C A AND B CONTAIN THE ADJACENCY MATRICES OF THE TWO GRAPHS, 
C C AND SC CONTAIN SPECTRA DUPLICATION INFORMATION. 

C P CONTAINS A PERMUTATION OF THE SET OF ORDER N-l SPECTRA OF A INTO 
C THOSE OF B. 

C MACHP IS THE PRECISION WITHIN WHICH EQUALITY IS TO BE TESTED. 
C MAXT IS THE MAXIMUM NUMBER OF PERMUTATIONS THAT COULD BE TRIED IF 
C THE OPTIMAL EIGENVECTOR CORRESPONDED TO AN EIGENVALUE OF 
C MULTIPLICTY ONE. 
C N IS THE ORDER OF THE GRAPHS. 
C NM IS THE DIMENSION OF THE ARRAYS. 
C 

C ON RETURN: 

C A AND B ARE UNCHANGED. 
C C AND SC ARE DESTROYED. 

C P CONTAINS A PERMUTATION OF THE COMPONENTS OF AN EIGENVECTOR OF 
C A SUBGRAPH OF A INTO AN EIGENVECTOR OF A SUBGRAPH OF B. IF 
C ERR=0 THEN THE PERMUTATION RELATES A AND B. 

C ERR = IF A PERMUTATION OF A INTO B WAS FOUND. ERR= 1 IF ALL 
C POSSIBLE PERMUTATIONS FAILED TO RELATE A AND B. ERR=2 IF ALL 
C PERMUTATIONS CHECKED FAILED BUT NOT ALL POSSIBLE PERMUTATIONS 
C WERE TESTED DUE TO THE RESTRICTIONS IMPOSED BY MAXT. ERR=3 IF 
C ALL EIGENVECTORS CORRESPONDING TO EIGENVALUES OF MULTIPLICITY 
C ONE HAD THREE OR MORE IDENTICAL COMPONENTS. ERR IS NEGATIVE IF 
C AN ERROR CONDITION WAS RAISED IN THE EISPACK SOFTWARE. 
C MACHP, MAXT, N, AND NM ARE UNCHANGED. 
C 

C WORKING VARIABLES: XA , XB , DA , DB , EA , EB , FA , FB , S D, SE , G ,H 
C 

INTEGER NM,N, MAXT, ERR, I , J , K , L ,M , LT , LM , PN, JA , G C NM) , HC NM) ,P(NM) , NN 

REAL A(NM,NM) ,B(NM,NM) ,C(NM,NM) ,XACNM,NM) ,XBCNM,NM) ,DACNM) , DB(NM) , 
+ EA(NM) ,£B(NM) ,FACNM) , FB( NM) , SC ( NM) ,SD(NM) , SE(NM) , MACHP , MACHPD 
+ ,EIGVAL,LB,UB,RD,FN 
C 

C FN IS THE FUNCTION FOR X- 1 IN THE 1,X FORM OF THE ADJACENCY MATRIX. 
C 

FN=N/2 

NN=M-i 

PN=10 00 
C 

C SEARCH FOR AN OPTIMAL EIGENVECTOR AMONG THOSE THAT CORRESPOND TO AN 
C EIGENVALUE OF MULTIPLICITY ONE. 
C 

DO 9 I=1,N 
C 

C UNIQUELY MARK THE I TH NODE OF A AND THE NODES CONNECTED TO IT. 
C 



DO 10 J=1,N 

DAC J)=AC I , J) 

AC J, J)=AC J, J)+DAC J) 

AC I , J)=3*FN 
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ACJ,I )=3*FN 
10 CONTINUE 

AC I, I )=AC I , I )+FN 
C 

C COMPUTE THE EIGENVALUES AND EIGENVECTORS OF THIS SUBGRAPH. A IS MOT 

C DESTROYED. 

C 

CALL TRED2 C NM , N , A , DB , FB ,XA) 
CALL TQL2 C NM , N , DB , FB , XA , ERR) 
ERR=-ERR 

IF (ERR.NE.O) RETURN 
C 

C RESTORE A TO ITS ORIGINAL VALUE. 
C 

DO 20 J=1,N 

ACJ, J)=AC J, J)-DA< J) 
AC I ,J)=DACJ) 
ACJ, I )=DACJ) 
2 CONTINUE 
C 

C SEARCH FOR THE OPTIMAL EIGENVECTOR OF THOSE JUST COMPUTED. 
C 

DO 80 J=l ,N 
LT=SCC I) 

MACHPD=MACHP*C 1 +ABSCDBC J) ) ) 
C 

C CHECK THAT THE ASSOCIATED EIGENVALUE HAS MULTIPLICITY ONE. 
C 

IF CJ.EQ. 1 ) GO TO 30 

IF CABSCDBCJ)-DECJ-1 )) .LT.MACHPD) GO TO 80 
IF CJ.EQ.N) GO TO 40 
30 IF CABSCDBCJ)-DB(J+1 )). LT.MACHPD) GO TO 80 

C 

C FIND AND RECORD THE NUMBER OF DUPLICATED PAIRS OF COMPONENTS. IF A 
C TRIPLICATE IS FOUND, DISCONTINUE CONSIDERATION OF THAT EIGENVECTOR. 
C 

40 DO 60 L=l ,NN 

M=L+1 
LM=Q 
DO 5 K=M,N 

IF (ABSCXACL, J)-XACK, J) ) .GT.MACHP) GO TO 50 
IF CLM.EQ. 1 ) GO TO 80 
LM=1 
LT=LT+1 
5 CONTINUE 

60 CONTINUE 

IF (LT.GE.PN) GO TO 80 
C 

C UPDATE CURRENT OPTIMAL EIGENVECTOR INFORMATION. 
C 

PN = LT 
JA=I 

EIGVAL=DBC J) 
DO 7 K= 1 , N 

SDCK)=XACK, J) 
70 CONTINUE 

C 
C IF THE EIGENVECTOR IS SUITABLE -NOT NECESSARILY OPTIMAL- TERMINATE 
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C THE SEARCH. 
C 

IF CLT.LT. 10 ) GO TO 100 

8 CONTINUE 

9 CONTINUE 
C 

C IF NO SUITABLE EIGENVECTOR COULD BE FOUND, TERMINATE THIS PORTION 

C OF THE ALGORITHM WITH APPROPRIATE MESSAGE. 

C 

IF (PN.LT. 1000) GO TO 100 
ERR = 3 
RETURN 
100 MACHPD=MACHP*C 1 +ABS (EI GVAL ) ) 
LB=EIGVAL-MACHPD 
UB=EIGVAL+MACHPD 
C 

C STORE THE SELECTED EIGENVECTOR IN THE FIRST COLUMN OF XA. 
C 

DO 110 I = 1,N 

XACI, 1)=SD(I) 
110 CONTINUE 
C 

C SELECTIVELY STORE THE PERTINENT INFORMATION HELD IN C AND SC IN ORDER 
C TO FREE STORAGE FOR FUTURE USE. 
C 

M=SCC JA) 

DO 120 1=1, M 

GCI)=CCJA, I) 
120 CONTINUE 
C 

C FIND AND EXAMINE THE CORRESPONDING EIGENVECTOR IN THE SUBGRAPHS OF B 
C WHICH HAVE SPECTRA MATCHING THE SPECTRA GENERATED FROM THE SUBGRAPH 
C OF A THAT PRODUCED THE OPTIMAL EIGENVECTOR. 
C 

DO 170 L=1,M 
C 

C COPY B INTO C. 
C 

DO 140 1=1 ,N 

DO 130 J=1,N 

CCI,J)=BCI ,J) 
130 CONTINUE 
140 CONTINUE 
C 

C FIND AND UNIQUELY MARK THE APPROPRIATE NODE IN B. 
C 

I=GCL) 

DO 150 J=l ,N 

CCJ,J)=C(J,J)+C( I, J) 
CCI,J)=3*FN 
CC J, I )=3*FN 
15 CONTINUE 

CCI,I)=CCI,I )+FN 
C 

C FIND THE DESIRED EIGENVECTOR - C I S DESTROYED. 
C 

CALL TRED1 ( NM , N , C , DA , DB , SC ) 

CALL BISECT C N , , DA , DB, SC , LB ,UB , 1 ,PN , RD , LM, ERR, SD, S£) 
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ERR«2*ERR-PN+1 

IF (ERR.NE.O) RETURN 

CALL TINVIT CNM t N, DA ,D5, SC , PN , RD , LM, FB , ERR, SD , SE , EA, EB , FA) 

£RR=-ERE 

IF (ERR.NE.G) RETURN 
CALL TRBAK1 ( NM , N , C , BB, PN , FB ) 
C 

C STORE THE EIGENVECTOR IN THE FIRST COLUMN OF XB. 
C 

DO 16 u = i , N 

XE(J, 1)=FB(J> 
16 CONTINUE 
C 

C SINCE BOTH EIGENVECTORS ARE STORED IN COLUMN ONE OF MATRIX ARRAYS, 
C SUBROUTINE SINGLE CAM EE USED TO GENERATE AND TEST UP TO MAXT 
C PERMUTATIONS, 
C 

CALL SINGLE ( A , 3 ,XA ,XB, DA , DB ,H ,P , ERR , MACHP , MAXT , i,N,NM) 
IF (ERR.EQ.O) RETURN 
IF (ERR.EQ.4) NN=N 
170 CONTINUE 
ERR=1 

IF CNN.EQ.N) ERR-2 
C 

RETURN- 
END 
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4. Computational Experience 

The UNI VAC 1108 computer at the National Bureau of Standards was used to test the computer 
program on well over 200 pairs of graphs specifically selected for use because of their difficulty. The 
results have been encouraging. A Brookhaven National Laboratories group has compiled [2] a list of 
149 pairs of nonisomorphic graphs of order varying from 7 to 10 that have 0,1 adjacency matrices with 
identical spectra. These cospectral graphs arise in connection with the energy field of a Heisenberg 
model ferromagnet. The computer program previously listed distinguished each of the 149 pairs in a 
total time of 3.9 cpu seconds or, equivalently, each pair of graphs was determined to be nonisomorphic 
in approximately 0.026 cpu seconds. Of the graphs, 109 failed the valence test, 28 the node duplication 
test, and the remaining 12 pairs failed to have 1,# modified adjacency matrices with identical spectra. 

With some modification the algorithm can be programmed to sort a large set of graphs into subsets 
characterized by node valences, node duplication, l,x modified adjacency matrix spectra, and in more 
difficult cases, submatrix spectra. Once this is accomplished, the groups in each subset can be pairwise 
checked for isomorphism using the program listed in this paper. This procedure was employed on a 
collection of 35 graphs ranging in order from 5 to 9. Many of the graphs selected were known to be 
isomorphic and all had large automorphism groups, thus making it as difficult as possible for the 
algorithm to find permutations or to determine that none exist. By means of a simplified form of the 
sorting procedure, the 35 graphs were divided into 19 subsets in approximately 5 cpu seconds. In this 
particular case, the spectra test was a sufficient condition for graph isomorphism. To be sure, however, 
permutations were generated between each pair of graphs in each subset. A total of 16 permutations 
were generated, 13 by means of SUBROUTINE SINGLE and 3 requiring SUBROUTINE DOUBLE, 
in a total time of 1.6 cpu seconds. 

None of the graphs yet tested, however, required the use of SUBROUTINE SUBVAL or SUB- 
ROUTINE SUBVEC. In response to a July, 1975 query about the existence of more challenging graphs, 
D. Corneil provided 16 pairwise nonisomorphic, order 25 graphs, any two of which required approxi- 
mately 60 seconds on an IBM 370 (no model number was mentioned) to be distinguished by then 
current methods. The graphs were 3-strongly regular and stochastic with node valence 12. None of the 
graphs has duplicate nodes and each had identical l,x adjacency matrix spectra. The submatrix spectra 
test, however, successfully sorted the group into 8 subsets, each containing two graphs. It is interesting 
to note that the two graphs contained in each subset were duals of each other. This initial sorting was 
accomplished in just 86 cpu seconds. When the algorithm was run on the 8 pairs of graphs, 5 were 
separated by means of the submatrix eigenvector test. The test was inconclusive on the remaining 3 
pairs of graphs. The total run time was 150 cpu seconds. Of this total, 86 cpu seconds was spent dis- 
tinguishing the 5 pairs and 64 cpu seconds on the unsuccessful effort to resolve the final 3 pairs. In 
addition, the algorithm was applied to 5 of the graphs and constructed permutations. The program pro- 
duced a permutation relating each of the 5 pairs of graphs by means of SUBROUTINE SUBVEC in 
81 cpu seconds. 

Approximately 90 percent of the time required by the program is used in the calculation of eigen- 
values and eigenvectors. Thus selection of optimal eigenvalue and eigenvector computing software is 
critical when time restrictions are sizeable. Since the EISPACK software referenced in the listing of the 
algorithm proved to be satisfactory for our needs, no significant effort was made to find more efficient 
software. 

Tests on large numbers of random graphs were not conducted as it was felt that such tests would 
yield little meaningful data. Experience gained from working with the eigenvectors and eigenvalues of 
many graphs showed that only a relatively minute number of pairs of graphs exist for which the 
question of isomorphism is at all challenging for the algorithm and no such graphs are likely to occur 
in even a very large sample of randomly generated graphs. For the average pair of isomorphic graphs 
of order n, it is estimated that a permutation will be generated in well under 10' 4 n B cpu seconds UNI- 
VAC 1108 time. Even faster times are achievable for random nonisomorphic pairs of graphs or through 
modifications in programming which first sort graphs by node valences, node duplicate structure l 9 x 
modified adjacency matrix spectra, and in special cases (such as Corneil's collection of order 25 graphs) 
submatrix spectra, and then use the algorithm described in this paper to further reduce the subsets into 
classes for which the graphs are pairwise isomorphic. 

Thus, with efficient application of the algorithm it is possible to divide a very large number of 
graphs into isomorphism equivalence classes in a reasonable amount of time. It is further possible to 
test pairs of very large graphs for isomorphism and if one exists, to generate permutations between two 
very large graphs. 
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